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Abstract 

We generalize the Kuramoto model of globally coupled oscillators to multifrequency communities. 
A situation when mean frequencies of two subpopulations are close to resonance 2:1 is considered in 
detail. We derive uniformly rotating solutions describing synchronization inside communities and 
between them. Remarkably, cross-coupling between the frequency scales can promote synchrony 
even when ensembles are separately asynchronous. We also show that the transition to synchrony 
due to the cross-coupling is accompanied by a huge multiplicity of distinct synchronous solutions 
what is directly related to a multi-branch entrainment. On the other hand, for synchronous popu¬ 
lations, the cross-frequency coupling can destroy a phase-locking and lead to chaos of mean fields. 
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I. INTRODUCTION 


Models in the form of coupled oscillators are ubiquitous in various scientific fields, ranging 
from physics and chemistry |1] to biology [2], as well as in some interdisciplinary applica¬ 
tions |3]. In many cases dynamics of oscillatory ensemble can be successfully studied in the 
phase approximation BE]. When the coupling between the oscillators is relatively weak, 
one can neglect changes in the amplitude dynamics of natural limit cycles of the oscillators, 
and describe the system in terms of the phases only. This technique is known as phase 
reduction, and it represents, basically, one of the few rigorous mathematical approaches to 
study complex non-equilibrium nonlinear oscillatory dynamics. 

The simplest setup here represents a globally coupled ensemble with weak interaction 
and relatively close natural frequencies. The phase reduction here leads to the system 
of globally coupled phase equations where interaction between the oscillators is described 
by 2tt - periodic function of phase differences H EHH]. The classical and well-studied 
Kuramoto-Sakaguchi model appears when one consider only the hrst Fourier mode in the 
interaction function what leads to simple sinusoidal coupling. There is almost 40 years of 
intensive studies dedicated to explanation of bifurcations and dynamics in this model [9] . A 
surprising recent result discovered a possibility of low-dimensional description of the classical 
Kuramoto model in terms of macroscopic order parameters [TU1 - 1T2] . However, this reduction 
to low-dimensional systems does not imply simplicity of dynamical behavior. In opposite, the 
authors ra report on quite complicated phase transitions and bifurcations in the Kuramoto- 
Sakaguchi models. 

The cases of multi-harmonic coupling functions mm appear to be more complicated and 
usually responsible for new dynamical effects in comparison to classical setup with purely 
sinusoidal function. In large ensembles the multi-harmonic case leads to appearance of so- 
called multi-branch entrainment modes with a huge multiplicity of possible synchronous 
solutions [S1EI[I5]- The latter also leads to non-trivial noise-induced effects |16j . 

One of the directions in this growing theoretical held is dedicated to multi-frequency os¬ 
cillator communities. As it was mentioned before, the Kuramoto-type models were obtained 
under assumptions of weak coupling limit and closeness of natural oscillator frequencies. 
However, when the distribution of the frequencies is huge in comparison to the interaction 
strength, the phase reduction leads to another types of phase models HMH]. A natural 
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setup here implies existence of a certain number of oscillator subpopulations (communities), 
such that the frequencies inside each population are close, but differ signihcantly across the 
distinct communities. This situation is inspired by theoretical and experimental results from 
neuroscience I2D1. indicating that distinct interacting brain areas exhibit different natural 
oscillatory rhythms. 

In this paper we consider a particular problem when distinct oscillatory communities have 
natural frequencies close to a high-order resonance. First, we derive general phase equations 
for globally interacting ensembles and distinguish different types of resonant coupling which 
may appear in the system. Next, we concentrate on the simplest case of two interacting 
population whose mean frequencies are close to an 2:1 resonance. The aim of the paper is 
to demonstrate on this simplest example, what one can expect from the effects high-order 
resonances. To describe the dynamics, we adopt the self-consistent approach developed 
in na for calculation of stationary order parameters for multi-harmonic coupling functions. 
Our analysis will show that the model exhibits reach dynamical behavior including multi¬ 
branch entrainment (multiplicity) and chaotic collective oscillations. 

II. PHASE EQUATIONS FOR RESONANTLY COUPLED POPULATIONS 

In this section we will present a general scheme of coupling in resonant, multifrequency 
populations of oscillators. We will assume that each oscillator is described solely by its phase 
0, which satishes the following equation 

^ = u + S{(j))F 

where u is oscillator’s natural frequency, S'(0) is its phase response curve, and F is the force 
acting from other oscillators. In order to simplify notations, we from the beginning will 
consider a thermodynamic limit, where the number of units in all populations and subpop¬ 
ulations tends to inhnity (although at the end we will also write the governing equations for 
a hnite size case). We assume that the ensemble is divided into M distinct subpopulations 
(we will use index n for referring to them), around M distinct mean frequencies uJn- Addi¬ 
tionally, there can be a small deviation from the mean frequency A (typically described by 
a unimodal distribution around zero). We now introduce slow phases, by writing explicitly 
fast rotating terms ~ a7„t. In fact, we can also chose frequencies of fast rotations to be 
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close, but not exactly equal, to ojn- We will use this freedom to be able to make perfect 
averaging below. Our slow phases = 4>n — ^nt satisfy equations 

0n('^) = A + Sn{^nt + ^n)Fn (1) 


where now individual mismatches A for the group n are distributed generally asymmetrically, 
with some small shift UJn — ^n- 

Next, we assume that coupling between the groups and inside each group is due to mean 
helds only. These mean helds for each subpopulation are represented by generalized order 
parameters 

where averaging is over the distribution of the slow phases following from ([^ and over the 
distribution of A. The introduced order parameters Z are slow functions of time as they 
are dehned via the slow phases: 

^^(n) ^ ('2) 

In general, the force acting on the oscillators of the group n is from all other groups, and 
is a nonlinear function of order parameters, which one can expand in powers of them. We, 
however, in this paper will restrict ourselves to the linear coupling only, i.e. we will assume 
that Fn is a linear function of order parameters: 


F ( yd) 7(2) 'i _ 


k,m 


,{m) y{m) 
^n,k ^k 






k,m 


Representing the phase response function Sn as a Fourier series 


( 3 ) 


p 

and substituting this in Eq. Q, we obtain 




'5n(A) — A + 


e JpvJp^nt 
onp^ c 




, /c,m 


p,k,m 


( 4 ) 


Now one has to perform averaging of Eq. Q, to reveal evolution of the slow phase. 
The fast terms on the r.h.s. are those containing explicit time dependence with one of the 
frequencies or with a combination of them. Such a combination can be small, this is 
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exactly the case of a resonance that is of special interest for us. Here, we use the freedom 


in the choice of particular values of to make the resonance exact. This means that 


some combination of frequencies vanishes exactly. Performing averaging means just 


keeping these terms on the r.h.s. of Eq. Q, and neglecting all other containing explicit time 
dependence. 


Expansion Q can be treated in many setups of particular resonant conditions, we describe 


here some evident cases: 

• One population of oscillators. In this case only one frequency O exists. Here the 
only terms surviving the averaging are those with p + A; = 0, this leads to the Daido 
model [8]. 

• Two subpopulations. Here the main interest is in the resonance of two frequencies 
Hi,f22. The simplest case is just the second-harmonic resonance: H 2 = 2Hi. In this 
case only those cross-population coupling terms with p-|-2fc = 0 survive. Similarly, for 
high-order resonances like aVL 2 = bQi (with integer a, b) the terms with ap + bk = 0 
contribute. 

• More than two subpopulations. One can see from (|^, that in the case of linear 
coupling, there is no direct interaction involving more that two subpopulations. So 
the resulting coupling is a combination of terms stemming from pairwise resonances. 

We restrict ourself in this paper to the simplest case of two resonant subpopulations with 
H 2 = 2Hi. As described above, after averaging only terms where combinations ~ (Hi^ 2 —f2i^2) 
and ~ (f 22 —2f2i) appear, survive, for the interaction within one and between subpopulations, 
respectively: 



k k 


( 5 ) 



k k 

We now insert here the dehnition of the slow order parameters ([^ and obtain 

0i(A) = A + 


k 


k 


( 6 ) 
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where averaging is over variables with tilde. Now we can dehne effective coupling functions 
inside subpopulations /ii ,/22 and coupling functions across subpopulation /i 2 ,/ 2 i as 

k 

k 

Now we can formulate equations for hnite populations, replacing () by corresponding sums. 
We assume that subpopulations 1 and 2 have Ni and N 2 units, respectively. Furthermore, 
one can now also transform back to the original fast phases, because in the averaged formu¬ 
lation the absolute values of the frequencies do not play any role. Denoting the phases in 
the subpopulation at a smaller frequency (we will also call it the hrst subpopulation below) 
as (j)p, and the phases in the subpopulation at a larger frequency (referred hereafter as the 
second subpopulation) as V’p, we get 

^ A'l ^ N2 

^ k=l ^ p=l 

^ N2 ^ Ni 

'^<1 = JfYl / 2 i ( 20 p - 'Ipq) 

^ k=l ^ p=l 

where we also have split notations for frequencies in two subpopulations. This system is 
a generalization of the Daido model [ 8 ] to two resonantly coupled ensembles. Below we 
will consider the case where coupling functions / contain the first harmonics only; this will 
correspond to the Kuramoto-Sakaguchi-type coupling. In this case each coupling function is 
determined by two parameters, the amplitude and the phase shift. One of the phase shifts 
in the cross-coupling can be set to zero by shifting all the phases in one subpopulation with 
respect to another one. Thus, our coupling functions will be: 

/ii(x) = £isin(x-ai), f 22 {x) = e 2 sm{x-a 2 ), /i 2 (a;) = 71 sin(x-/?), /22 (a;) = 72 sin x 

Next, we £x the distributions of the frequencies. As after the averaging the system 
is invariant under transformation 0 —)■ 0 -|- 2At, 'ip —)■ p) + At for arbitrary A, we can 
set the average value of the natural frequencies in the hrst subpopulation 0 to zero, the 
average frequency 6 in the second subpopulation is the relevant parameter responsible for 
the mismatch. We will assume the frequencies to be distributed according to Lorentzian 
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distributions, with equal widths. Because we still have a freedom of changing the time scale, 
we will assume that this width is one: 


9i{^) = 


92(1^) = 


7r{uj‘^ + 1) 7 r((z/ — 5)2 + 1 

The resulting microscopic system of oscillators to be considered below reads 


(9) 


£1 


Ni 


= + sin( 0 fc - - ai) + ^ “ Z^) 

^ k=l 2 


N 2 


£2 


N 2 


^ sin(V’;t 

iVo 
^ k = l 


A 


" 2 ) + 

^ k=i 


Ni 


( 10 ) 


t/’m) 


( 11 ) 


with frequencies dehned according to the distributions 

We now also write down the basic equations in the thermodynamic limit. Here three 
complex order parameters Xi,X 2 , Y appear dehned as 

Xfc = = jj d^dw k = 1,2 , 

Y = = Jj dijdu g2{v)p{'ip\A£"^ 

while equations for the phases are 

0 = a; + EiXi sin( 0 i — 0 — ai) + 71 Y sin( 0 j^ — 20 i — (d) 

-tjj = ly + e 2 Y sin( 0 j^ - -0 - 0 : 2 ) + 72-^2 sin (02 - A 
The formulated system of equation will be subject of our analysis below, where we will 
concentrate on main dynamical effects caused by resonant cross-coupling. In numerical 


( 12 ) 


simulations we will use microscopic equations ( 10 ), while in the theoretical construction the 


thermodynamic limit formulation (ITp!2) will be used 


III. SELF-CONSISTENT SOLUTIONS IN THE THERMODYNAMIC LIMIT 


Here we will present the self-consistent scheme allowing us to hnd stationary (or, more 
generally, uniformly rotating) synchronous solutions of the system (II][l2). 


A. Ott-Antonsen ansatz for the second subpopulation 

The problem partially simplihes by the observation, that for the subpopulation at 
the double frequency the Ott-Antonsen ansatz m can be applied. Indeed, the second of 
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eqs. (12) can be rewritten as 


^/^ = i/ + Ini(H(t)e-'^), H = £ 2 e-'“^Y + 72 X 2 (13) 

According to the Ott-Antonsen theory, equation for the order parameter Y obeys (un¬ 
der some additional assumptions which we assume to be satished here), in the case of a 
Lorentzian distribution (|^, an ODE 

Y = Y(2(5-1)-^(H*Y2-H) = Y(i(5-l)-^(e-*“^|Y|2-e*“^)-y(X;Y2-X2) (14) 


B. Uniformly rotating ansatz 


We now construct solution for the ensemble of oscillators ip. Here we cannot use the 
Ott-Antonsen ansatz, because the latter is only applicable for the driving terms possessing 


one harmonics of the phase, like in (13). The equations for ip possesses both the first and the 


second harmonics. In order to find stationary values of the mean fields, we will adapt the 
self-consistent scheme developed in Refs. [15] for the deterministic bi-harmonic Kuramoto 
model (for the noisy case a similar method can be used, see [161 El])- 

In this self-consistent approach one finds uniformly rotating distributions, i.e. distribu¬ 
tions that are stationary in a rotating reference frame. Let us denote the frequency of this 
frame D, it will be determined self-consistently as a result of the calculations. According to 
this, we introduce constant phases of the order parameters 


0^ = Uf, 02 = 02 - 2Dt, 0=0- 2VLt 


( 15 ) 


(here the phase shift of the hrst order parameter Xi is set to zero, this can be always done 
by the time shift). Also, we introduce a new phase variable 

ip = (j) — VLt + ai 


distribution of which is expected to be stationary. This variable obeys 


ip = oj — VL + SiXi sin(—(^ 9 ) -|- 7 iY sin(0y -|- 2 q!i — j3 — 2ip) 


( 16 ) 





C. Stationary solution in a parametric form 


To proceed with self-consistent solntion, it is convenient to introdnce fonr auxiliary pa¬ 
rameters {R, u, V, z} = 'P in the following way: 


SiXi = i?sinM, 7 iy = Rcosu, hi = zR^ v = 9y + 2ai — {3 


( 17 ) 


Now (16) takes the following form: 


(p = R{x — z — sinM sin (p — cos nsin( 29 ? — n)) = R{x — z — h{u, v, ip)) 


( 18 ) 


We denoted x = u/R and h{u,v,p) = sin n sin (p -1- cosMsin(2(p — v). At some constant 


values of parameters P in (18), at each value of x one can hnd stationary distribution 


function p((p|x,P), and then calculate the corresponding complex order parameters: 


Xi = e-^^^R Jj p{p\x,P)R^g{Rx)dxdp = 

jj p{p\x, Py"^^g{Rx)dxdp = RF 2 {Py^^^'^'> (19) 

= [[ dxdpp{p\x, Py^^g{Rx) , m = 1, 2 . 


Our next goal is to calculate the integrals F^iP), for this we need to hnd, using the 


dynamical equation (18), the stationary distribution function p{p\x, P). Let Hmin and Hmax 
denote the global minimum and the global maximum of function h{u,v,p), correspondingly 
(FigQb)). All the oscillators can be separated into locked ones (for Hmax > |a) —> Hmin) 
or rotating, unlocked ones [x — z > Hmax oi x — z < Hmin)- The distribution function of 
rotating oscillators (index r) is inversely proportional to their phase velocity: 

C{x) 


Pr{p\x,P) = 


\x — z — h((p,M,n)| ’ 


( 20 ) 


where C{x) is the normalization constant: 


C{x) = 


tStt dtp 

Jo 


/O \x-z-y\ 

The stationary phases of locked oscillators (index 1) can be found from the following 
relation: 

X — z = h{u, V, p) . (21) 
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FIG. 1. (a) Regions Vi and V 2 in the plane of parameters {u,v): Domain Vi corresponds to a 
double-well form of function h{u,v,(p) (Fig. [^b,d)), while in V 2 h{u,v,(p) has a single-well form 
like sown in Fig. [^c). (b) Example of function h{u,v,ip) with 4 extrema is presented. There 
are two stable branches (solid curves) for stationary phases of locked oscillators. The left branch 
If = 'I'i(a:,P) is larger than the right one ip = 'I' 2 (a:,P). {ipi, 2 ,xip) denote coordinates of the 
extrema corresponding to the branch 'I'l, while ((/ 53 , 4 ,^ 3 , 4 ) denotes extrema at \k 2 - In the domain 
h{(p) G there is a bistability on the microscopic level: in this domain the oscillators can 

be locked either on the branch in the range ip G [ip\,ip 2 ] or on the branch 'I '2 in the range 
ip G [v:’3,<^4]. (c) Example of function h{u,v,(p) with only two extrema and one stable branch 
ip = 'I'i(x,P) (solid curve). 


When finding ip as a. function of x for non-rotating (locked, index /) phases, we have to 
satisfy an additional stability condition > 0 that follows from the dynamical equa¬ 

tion (18). In the {u,v) plane there are two regions Vi and V 2 (Fig. [^a)) with qualitatively 
different properties of the system (18) and different types of distribution function pi{ip\x, P), 
corresp ondingly: 
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a. {u, u} G Vi In this case function h{u, v, ip) has a double-well form like shown in 
FigQb). According to ( [I^ , oscillators can be located on two possible stable branches 
highlighted by solid curves in FigQb): the first branch is p = P) in the range p G 

[PI)Pq\ and another branch is (p = \E' 2 ( 2 :, P) for p G [(^ 3 ,(^ 4 ]. Here and below we assume 
P) to be the biggest stable branch. In the range {x—z) G (xi, x^) (Fig.j^b)) there is an 
area of bistability on the microscopic level: the oscillators with the same natural frequency 
X can be locked at two different phases \I'i(x,P) and \I' 2 (a^, P)- Therefore, the distribution 
function has the following form: 


Pii^\x,P) 


' (1 - - Ti(a:, P)) + S{x)S{p - ^ 2 ( 0 :, P)) 

for {x — z) & [x\^ , 

5{p - ^ 1 ( 0 ;, P)) for {x - z) e [xi, X 2 \ \ [x\,X 2 ] , 
^ 6 {p- T 2 (x, P)) for {x- z) e [x 3 , X 4 ] \ [x\,xl] . 


( 22 ) 


Here 0 < S{x) < 1 is an indicator function describing the redistribution over the stable 
brunches; this function is arbitrary. 

b. {u,v} G V 2 In the second case, function hiu^v^ip) has only two extrema (Fig. fjc)) 
and there is only one stable branch p = \['i(a:, P). The distribution function is: 


pi{p\x, P) = 5{ip - \I'i(a;,P)) for x e {z + xi,z + X 2 ) ■ (23) 


Taking into account the obtained expressions for the distribution function ( [^22[23[ ), the 
integrals in (19) can be rewritten as a sum of five terms: 


F^(P)e*«-(P) = / dpe^^^g {R{z + h)) - - 

dp 

r^4: 




dh 


dh 




dpR^^S{z + h)g{R{z + h)) — + / dpR^^g{R{z + h))- - 


dp 


-v'i 


' ¥>3 

dh 




dpR^^ (1 - S{z + h)) g {R{z + h))^ + 



2tt 


dxdp 


X Jo 


dp 

g{Rx)C{x)R^'^ 
\x — z — h\ 


(24) 


Here the first and the second terms stand for integration over the first branch \ki in the range 
[Pi,Pq\. The second term accounts for certain redistribution S{x) of oscillators between 
the branches in the range [p\-,p^ 2 \ (Fig- [^F)). Similarly, the third and the fourth terms 
correspond to integration over the possible stable branch '^2 in the range [(^ 3 ,(^ 4 ]. In the 
same way, the fourth term accounts for redistribution of oscillators between branches in the 
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range [ 99 I, (^ 4 ] (Fig. [^b)). In the last term the interval X = (— 00 , z + Hmin) U('^ + Hmax, 00 ) 
is the domain of frequencies where the oscillators are not locked. 


Now, using the integrals (24), one can calculate the absolute values of the complex order 


parameters Xi ^2 and the frequency hi as functions of introduced auxiliary parameters R, u, 


Xi, 2 (P) = i?Pi, 2 (P), 02 = Q 2 (P) - 2ai, fi(P) = Rz . 

Also, from the relations (iTpQ) one can conclude that the following holds: 


sinti ^ Rcosu ^ ^ 

/x-.\ ) = Ql(P); 7l = - ^7 -) P = 6y + 2Ql(P) — V. 


FAP) 


Y 


(25) 


(26) 


D. Accounting for coupling between subpopulations 


As one can see from the latter relations, the parameters of the internal interaction inside 
the hrst community £1 and ai are determined only by the parameters P. However, the 
constants of the cross-coupling 71 and {3 require calculation of the order parameter Y. 


Taking into account the transformation of variables (15), the uniformly rotating solution of 


the Ott-Antonsen equation (14) for the second population, the mean held Y is determined 


according to the following relation: 


Ye*®"(i(2H-h) + l) + 


e2Ye 


iOu 




72-A2 , 


*02 y 




= 0 . 


(27) 


This complex equation determines Y and Qy as functions of all other parameters, substitution 


of these values to Eqs. (26) will give the values of cross-coupling parameters 71 and /I. 


In general case solution of (27) can not be represented in an analytic form and one should 


use certain numerical methods to hnd them (a parametric representation of solutions may be 
possible, but we already have four auxiliary parameters, introducing another two appears 


not practical). However, in two special cases equation (27) can be reduced to a simple 


polynomial equation with analytic solutions available. Namely, (i) for £2 = 0, the problem 
reduces to a complex quadratic equation, and (ii) for the special case H = 5 = 0 and n = 0 , 


equation (27) reduces to a real cubic equation. The latter case corresponds to the simplest 


situation when there are no phase shifts in coupling functions a \^2 = /3 = 0 . 

Summarizing, the self-consistent approach for calculation of stationary synchronous solu¬ 
tions of the problem (ITp)2) consists of the following steps: (i) for a given set of parameters 


12 















P and indicator function S{x), one constructs the distribution function p{(p\x) using mi¬ 
croscopic dynamics ( 18 ). (ii) Next, using the function p{p\x) and equations ( ^ 25 [ 26 ), one 
determines the stationary values of order parameters Xi 2, rotating frequency and corre¬ 
sponding coupling constants ei and ai. (iii) In the following step, one should solve equation 


( 27 ) for any hxed values of £2, 0^2, 72 and 6. As a result, one get the stationary value for the 


mean held Y and remaining constants of the cross-coupling 71 and /3 from ( 26 ). 

The solution is in the parametric form: varying the set of auxiliary parameters P, together 
with £2, 0^2 and72, one gets different solutions for the mean helds Xi_2, Y, together with 
their dependence on the coupling constants £1, 71 and ( 3 . This can be done for any indicator 
function S{x), which determines re-distribution of the phases of the hrst subpopulation 
between possible stable locked states, if multi-branch entrainment is possible. 

In the next sections we will apply the self-consistent scheme to characterize main types 
of synchronous states existing in the system of two coupled subpopulations ( IIp!2 ). We will 
focus on the effects caused by the resonant cross coupling between the population, therefore, 
for the internal coupling we will consider the simplest situation when 


ai,2 = 0. 

For the sake of simplicity, we restricts ourselves to the following parameters area: 

ei= 82 = e and 71 = 72 = 7. 

It appears that the latter choice of parameters simplify the presentation of the results, 
nevertheless, it contains all the main effects peculiar for the high-order resonant interaction. 


IV. INTERNALLY ASYNCHRONOUS POPULATIONS, APPEARANCE OF SYN¬ 
CHRONY DUE TO RESONANT COUPLING 

We start with the analysis of the case when the populations are internally asynchronous, 
hence, without cross-coupling 7 = 0 the only stable state for each population is asynchrony 
when all mean helds vanish Xi^2 = 0 , Y = 0 . For the Lorentzian distribution of frequencies, 
the synchronization sets in at the critical coupling e = 2 . Therefore, in the following section 
we will concentrate on the case e < 2, i.e. the internal coupling inside each population is 
insufficient to maintain synchrony in the system (or even repulsive). The frequency mismatch 
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(c) 



(d) 



FIG. 2. (a) The surface M depicts the boarder of synchronous states in the parameter space 
(h,/ 3 , 7 ): above the surface synchrony with Xi ^2 7 ^ 0 and T 7 ^ 0 exists, below only asynchronous 
state is possible. The internal coupling e = 1 for each population, hence, population are internally 
asynchronous, (b) The dependence of order parameters on the coupling constant 7 is shown for 
6 = (3 = 0 and e = 1. The curves denote theoretical calculations using self-consistent scheme, 


markers correspond to the direct numerical calculations of the finite-size ensemble (10) for N = 


8 X 10^. (c,d) Cuts of the surface M are shown for constant values of the frequency mismatch 6 
(in the panel (c)) and the phase shift /3 (in the panel (d)). 


6 together with cross-coupling constants 7 and phase shift /3 constitute a set of main control 
parameters in the system. 

Figure a) shows the area of existence of stationary synchronous solutions in the 3- 
d parameter space {5,(3,^). The surface M depicted in the Fig. [^a) denotes the border 
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of existence of synchronous states: above the surface there exist stationary synchronous 
solutions with Xi 2 7^ 0 and V ^ 0, below M only asynchronous state exists and is stable. 
Figurej^b) explains the bifurcation diagram depicted in Fig. j^a), here we £x the parameters 
S = /3 = 0 and plot order parameters Xi ^2 and y as a function of parameter 7 (the latter 
corresponds to the vertical line passing through the origin in the Fig. |^a)). As one can see 
from the plot, there is a minimal critical coupling 'fcr corresponding to the point m where 
two branches of synchronous solutions arise. The upper branch appears to be stable, what 
is confirmed by direct numerical simulation of the finite-size ensemble. The lower branch 
is unstable and disappears in the point z, merging with the trivial state. Apparently, the 
family of the points m obtained at different values of and 7 constitutes the surface M 
depicted in the Fig. |^a). 

The form of the surface is invariant under transformation 6 —)■ —S, (3 —)■ —/3, that is why 
only the part with h > 0 is shown in the Fig. |^a). Expectedly, M has a global minimum 
at the point 6 = f3 = 0 what means that, substantially, phase shift and frequency mismatch 
act against synchronization. Figures [^c,d) show several cuts of the surface M at constant 
values of (5 = const (in the panel (c)) and f3 = const (in the panel (d)). For the most part of 
the parameter range, the phase shift acts against synchronization, as one can easily see from 
the Fig. [^c ) where the borders of stationary synchronous states are plotted on the plane 
(/3,7). When the frequency mismatch mismatch is absent (h = 0), the curves are symmetric 
with respect to the line /d = 0 and critical coupling increases with growth of the absolute 
value of (3. However, it is not always a case for non-zero frequency mismatch. The examples 
for h 7^ 0 in the Fig. [^c ) clearly indicate a nontrivial fact that the global minima of the 
curves in the {(3, 7) plane are shifted towards negative values of f3. Similarly, on the (5,7) 
the boarder of synchronous states has global minima at non-zero value of S for finite phase 
shift (3 = TT/A. 

Remarkably, the transition to synchrony here is always accompanied by the multiplicity of 
different synchronous states with multi-branch entrainment [HIII3 in the hrst subpopulation. 
The issue of multiplicity for the bi-harmonic Kuramoto model was studied in detail in [T5] . 
The multiple synchronous states appear as a result of strong second harmonic ~ 6*2*^ in 
a global force acting on oscillators of the hrst subpopulation. Apparently, in order to get 
a synchronization in the ensembles due to the cross-coupling, the constant 7 has to be 
strong enough, as one can easily see from the bifurcation diagram in Fig. [^a). The latter 
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FIG. 3. The areas of existence of stationary synchronous solutions on the parameter planes (/3,7) 
are shown for the case e = 1 and <1 = 0. Different curves correspond to boarder of synchronous 
states with different indicator functions S{x) = a = const (different multi-branch entrainments). 
Above the curves solutions with Xi ^2 7^ 0, T 7^ 0 exist. Inset shows boundaries of synchronous 
states plotted for different values of constant e. From bottom to top e = 1.5, e = 1.0, e = 0.5. (b) 
Dependences of order parameters Xi^ 2 ^ Y on cross-coupling 7 are shown for states with different 
a (see legend). Solid curves denote solutions of self-consistent equations, markers denote direct 
calculations of the finite-size ensemble. Other parameters are: e = 0.5, /3 = 0, J = 0. (c) The same 
as (b) but for e = 1.5. 


implies that the coupling function h{u,v,ip) (see eq. (18)) always has a double-well form, 
hence, there is always a possibility to redistribute oscillators between two stable branches 
in different ways (in other words, to choose an arbitrary indicator function S{x) in the self- 
consistent scheme). As a result, for the case of internally asynchronous populations (when 
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e is not large enough), a family of synchronous states appears with distinct multi-branch 
entrainments. Figure |^a) shows the critical couplings 7 when synchronous states with 
distinct redistributions S{x) appear. For the sake of simplicity we chose S{x) = a = const. 
The dependences of order parameter Xi^2, F' on the cross-coupling 7 for different types of 
multi-branch entrainments (characterized by constant a) are presented in the Fig. [^b). The 
main state a = 0 arises hrst (i.e. at a minimal coupling strengths 7) in comparison to other 
states with multi-branch entrainment a 7^ 0 . Expectedly, in all cases increase of the coupling 
7 leads to increase of order parameters values, hence, more oscillators are entrained in both 
populations. 

Here in this section we have concentrated on the effects caused by the cross-coupling 7 
and paid less attention to the role of the internal coupling e. It is worth mentioning that, in 
the simplest form (pure sinusoidal coupling) the interaction inside the communities makes 
a relatively straightforward effect. Namely, increase of the coupling £ leads to enlargement 
of the area of synchrony existence in the parameter space (see inset in the Fig. [^. 


V. INTERNALLY SYNCHRONOUS POPULATIONS: CHAOTIC DYNAMICS. 


In this section we consider the case when e > 2 , hence, the populations are internally 
synchronous even without cross-coupling 7. Here we report on non-trivial phenomena when 
resonant cross-coupling introduces chaotic collective oscillations into the system. 

Figure |^a) shows the diagram of stationary synchronous states versus phase shift in 
the cross-coupling function f 3 . The solid curves correspond to solution obtained from the 
self-consistent approach above, while markers denote direct numerical calculations of the 


ensemble (10) at the same parameters. As one can easily see, the stationary states remain 


stable until [3 is less than certain critical value (indicated by colored area in the Fig. |^a)). 
However, when the phase shift becomes relatively close to tt, the synchronous solutions 
loses stability and immediately the system switches to a chaotic oscillation mode. The 
corresponding time series is presented in the Fig. |^b). Remarkably, chaotic oscillations 
are characterized by a drift of the phase difference ©2 — Qy (see the lowest panel in the 
Fig. gb)), so, the ensembles suddenly become unlocked from each other when entering to 
chaotic mode. 

Fig.g c) is aimed to explain the structure of parameter area where chaotic mode exists. 


17 




1.0 

0.0 


1.0 




wprpr 



1.0 

Yo.s 

0.0 


3 


(a) 


—TT+O.l —7r+0.2 —77+0.3 

(3 



(c) 




(d) 


7 


FIG. 4. (a) Dependence of order parameters on phase shift f3 is presented. Solid and dashed curves 
denote solutions obtained from the self-consistent approach. Stable solution corresponds to solid 
line, unstable to dashed line (stability was checked by direct simulation of the ensemble p^). 


Markers correspond to simulation of finite-size ensemble (10) for N = 10‘^ oscillators. The colored 


area denotes chaotic region with large amplitude of order parameters oscillations. Parameters 
e = 4.5, 7 = 2.8. (b) Time series of the hnite-size ensembles in the chaotic regime. Parameters: 
e = 4.5, 7 = 2.8, /3 = —3.0. (c) The region of existence of chaotic mode on the parameter plane 
(e,7) is presented for /3 = vr — 0.1, <5 = 0. (d) Dependence of the order parameter Xi is presented 
for e = 5, j3 = vr — 0.1 and <5 = 0. The gray area denote maximal amplitude of oscillations obtained 
from time series after certain transient period. 


As one can see, for each sufficiently strong internal coupling e, there is always a certain 
range of cross-coupling constant 7, where oscillations are irregular. With an increase of the 
cross-coupling 7, the system pass from area A (the area where regular synchronous solution 
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exists) to area B which corresponds to chaotic motion. As has been mentioned above, area 
B is characterized by a drift of the phase difference and large amplitude irregular oscillations 
of the order parameters (see Fig. |^d)). Further increase of constant 7 leads back to a regular 
stationary synchronous solution (Fig. |^d)). The size of the area B on the (£,7) plane is 
strongly related to the phase shift in coupling function: the closer the parameter /3 to vr, the 
larger the area B on the {e, 7) plane. 


VI. CONCLUSIONS 

The phase reduction is one of the few mathematical techniques which allows one to 
perform analytical studies of complex nonlinear oscillatory systems. Perhaps, the most 
popular and well-studied phase model is the classical Kuramoto system which describes 
ensemble of globally coupled oscillators with sinusoidal type of interaction function. The 
derivation of various Kuramoto-type models is based on the assumption of closeness of 
natural oscillatory frequencies. However, in many realistic situations oscillators may have 
definitely different frequencies; an example of this are neural populations that can produce 
brain waves wide across the spectrum. For multifrequency populations one has to extend 
typical model of phase equations; in previously considered cases such an extension also led 
to new dynamical regimes [TSl HH] . 

In the present paper we developed an extension of the phase synchronization theory for 
multifrequency resonant oscillator communities. After analysing general possible resnonant 
terms in linear mean-field coupling, we focused on the simplest high-order resonant case, 
when two communities of oscillators are globally coupled and have natural population fre¬ 
quencies close to the rational relation 2:1. First, given the assumption on mean population 
frequencies, we derived the simplest form of phase equations for high-order resonant interac¬ 
tion between two globally coupled communities of oscillators. Basically, the structure of the 
model consists of two main parts: the first part represents classical sinusoidal term describ¬ 
ing Kuramoto-type interaction inside each community; the second component has different 
form, it represents the resonant cross- coupling between the populations. Next, we combine 
two approaches described in Refs. na and in Refs. m to derive a self- consistent scheme 
for calculation of stationary synchronous solutions of the system in the thermodynamical 
limit. 
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In this paper we focused on investigation how cross-coupling promotes synchronization, 
and looked for novel dynamical effects due to the high-order resonance. Hence, we have 
considered two qualitatively different cases, in the first case the populations were internally 
asynchronous, so, the internal coupling strength was relatively weak. Here we constructed 
the bifurcation diagram showing how synchronous regimes appear in dependence on main 
parameters of the model. We demonstrated, that strong enough resonant cross-coupling 
results in stationary synchronous solutions appearing in both subpopulations. Thus, the 
synchrony can be only mutual. The nontrivial fact here is that the transition to synchrony 
due to the cross-coupling is always accompanied by multiplicity of distinct synchronous 
states, similar to the case of bi- harmonic Kuramoto model ng. In the second setup, we 
considered an opposite situation, when the internal coupling is strong, such that almost all 
oscillators are locked to the mean fields in the absence of the cross- coupling. Here we report 
on a quite non-trivial effect, that the cross-coupling can destroy the stationary synchronous 
state introducing chaos into the system. Mean fields of two subpopulations not only vary 
their amplitude chaotically, but the subpopulations also desynchronize from each other in 
the sense, that the phase shift between the mean fields is no more a constant, but performs 
a biased random walk. 
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